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Abstract: We apply the equivalent theory to orthorhombic anisotropic materials and provide 
{SJ ■ a general unit-cell design criterion for achieving a length-independent retrieval of the effective 

CN \ material parameters from a single layer of unit cells. We introduce a graphical retrieval method 

and phase unwrapping techniques. The graphical method utilizes the linear regression technique. 
Our method can reduce the uncertainty of experimental measurements and the ambiguity of 
phase unwrapping. Moreover, the graphical method can simultaneously determine the bulk 
values of the six effective material parameters, permittivity and permeability tensors, from a 
single layer of unit cells. 
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1. Introduction 

Metamaterials (MMs) are artificial materials engineered to achieve unusual electromagnetic (EM) properties that are 
not normally found in nature HI -|3). A comprehensive review can be found in the book|4|. Distinguished from photonic 
crystals, metamaterials are made from periodical structures with unit cells much smaller than the wavelength of light. 
A general method to construct large area 3D MMs is layer-by-layer fabrication technique|5]-(|8)- In stacked MMs, 
interaction between adjacent layers makes it difficult to extract bulk material parameters from a single unit-cell layer. 
It has been found that the retrieved effective metamaterial parameters are often dependent on the number of unit cells 
along the propagation direction[5, 9|. It appears that there is no clear methodology on how to accurately predict the bulk 
values of the effective permittivities and permeabilies through a single layer of unit cells. Currently there is no simple and 
effective way to resolve phase ambiguity in determining the phase of the transmission and reflection of electromagnetic 
fields because of phase wrapping. Although the phase ambiguity is a common issue for the parameter retrieval of the 
general composite materials, this issue becomes more significant for metamaterials where typically resonances, positive 
refractive index, and negative refractive index are all present in the same frequency band. In this paper, we apply 
Herpin's equivalent theorem |[T0l to orthorhombic anisotropic media and provide a simple way to accurately predict 
the effective bulk material parameters from a single layer of unit cells. We introduce a graphical retrieval method and 
phase unwrapping techniques, which can simultaneously determine the six material parameters, the permittivity and 
permeability tensors, from one unit cell. 



2. Equivalent theory 

Layer-by-layer fabrication method renders metamaterials intrinsically anisotropic. In MMs one unit-cell layer is usually 
composed of several sub-layers of different materials or nanostructures. We limit our discussion in nonchiral materials. 
To illustrate the key point we adopt multilayer method that has been used to model MMs by many groups lfTTI - lfT4l . We 
further assume that the principal axes of the sub-layers are parallel in each direction. In the principal coordinate system, 
the permittivity and permeability tensors are given by 















liny 





(1) 



where n = 1 , 2, • • • . The scalar terms e„j and jj, n j (j =x,y,z) are complex. Consider a monochromatic wave of frequency 
co with time dependence exp(— icot ) propagates inside the orthorhombic anisotropic materials. In each layer we have 

Vx(^ 1 -VxE) = k 2 (e n -E), 



Vx(£- 1 -VxH) = k 2 (fi„-H), 
where ko = co/c. If the plane of incidence is one of the crystal planes, the TE and TM polarizations are decoupled. 



(2) 




Fig. 1. A schematic shows how to construct a symmetric unit cell, which is composed of 
two original asymmetric unit cells. 



According to Herpin's equivalent theorem lfTol . every general multilayer is equivalent to a two-homogeneous-layer 
system and every symmetric multilayer is equivalent to a single homogeneous layer, characterized by an equivalent 



index and equivalent thickness fl5l . In metamaterials, a single unit-cell layer can be considered as a multilayer system. 
So it is equivalent to a two-homogeneous-layer system, denoted as AB. We can then construct a symmetric unit cell by 
cascading two unit cells as ABBA, which is equivalent to a single homogeneous layer according to Herpin's theorem. 
This process is illustrated in Fig. Q] Applying this methodology, a symmetric unit-cell layer can often be constructed 
regardless the number of sub-layers and complexity of each sub-layer in the original unit cell. The permittivities and 
permeabilities retrieved from the symmetric unit cell, which is composed of two original asymmetric unit cells, will 
represent the bulk material parameters. Thus, a length-independent description can be achieved. 

To make it more clear, let M represents the characteristic matrix of one symmetric unit cell. According to Herpin's 
equivalent theorem, it can be replaced by an equivalent single layer. Assume the x-z plane is the plane of incidence. For 
TM mode, i.e. H = (0,// v ,0) and E = (E x ,0,E z ), we have 

Mn =M 22 = cosv/ e , 

sm\j/e . (3) 

where y/e is the equivalent phase thickness of the symmetric unit cell; 2f e is the equivalent impedance. If the material 
contains AMayer symmetric unit cells, the characteristic matrix of the material is given by 

cosy* Jx sinVe \ =1 cos (^^) Tgia*(N\pe, . 
-/J^sini/Ze cosi/^ J \—i3? e sm{Ny e ) cos(Ny e ) 

A similar expression for TE mode can be obtained by replacing the 3? e with the negative admittance — W e . Here the J^ 
and W e are, respectively, the generalized impedance and admittance because they include incidence angle, i.e., Eq. <j4j is 
valid for both normal and oblique incidence. They are given by 

TM: 3- = _*_ kl= kle x Hy - -*5 , 

T ' " a ' (5) 

where (e x ,e y ,e z ) and (Hx,Hy,Hz) are, respectively, the effective permittivities and permeabilities of the equivalent layer. 
As shown in Eq. (@), the total phase of a AMayer system equals N times the phase \j/ e of the single layer; whereas 
the impedance 2f e or admittance W e is independent of the number of layers. These properties imply that the material 
parameters retrieved from the N layers of unit cells are the same as those retrieved from one layer of unit cells. In other 
words, the bulk metamaterial parameters can be predicted from a single symmetric unit cell. Note that the characteristic 
matrix of a N-layer asymmetric unit cells does not have the nice properties as those in Eq. (0J, As a consequence, the 
retrieved material parameters will be dependent on the number of unit cells along the propagation direction. Hence, the 
bulk material parameters cannot be predicted from one unit-cell layer. In other words, a length-independent description 
cannot be achieved for asymmetric unit cells. 

3. Graphical retrieval method 

In this section, we will introduce a graphical retrieval method based on Herpin's theorem and the previous method! 1 1 ]■ 
From above discussion a symmetric unit cell can be represented by three equivalent layers ABA. Let the z-direction 
perpendicular to the plane of the layers. Assume each layer is orthorhombic anisotropic with the principal axes parallel 
in each direction and the effective permittivity and permeability tensors are given by Eq. (HJ. Let the x-z plane be the 
plane of incidence. The characteristic matrix of TM mode is defined as 

M 2 \ M 22 





where 'i' refers to the input; and 'o' refers to the output. By matching boundary condition, the relationship between the 
scattering and characteristic matrices of a unit cell is given by 



M\ i = M 2 2 



{l+S 2l )(l-S l2 )+SnS 2 2 



M l2 



Mix = 



2S n 
l+S 2 i)(l+S n )-S n S 22 



(7) 



2S7T 



2Sn2f 
[l-S 2l )(l-S l2 )-S n S 22 



where 3fj and 3f are, respectively, the generalized input and output impedance whose values depend on the incidence 
angle and background permittivity and permeability at the input and output sides. In our definition, Stj(i = j) is related 
to the transmission and Sij(i ^ j) is related to the reflection. That means, for the TE polarization Sn and S 2 i are, 
respectively, the transmission and reflection coefficients of the electric field; whereas for the TM polarization the 



transmission and reflection coefficients of the electric field are given by 



£2^1; 



S\\ and —52i, respectively. Note that 



as long as the layer has inversion symmetry, M\\ = M 22 regardless the input and output impedances. For scattering 
matrix, S\\ = S 22 only if ^7 = 2f . In addition if M\i = M 22 , we then also have S12 = S21. From Eq. (01, the dispersion 
relation and the effective impedance 2f e of the unit cell ABA are given by 



cos(K e d) ~M\\= cos(ki z di) cos(k 2z d 2 ) — Tj + sin(^i z iii)sin(^2z^2) 
2 _ Mn = 2 sin(^ lz <ii ) cos(k 2z d 2 ) + r\ ^ 
e M l2 1 sm(ki z di)cos(k 2z d 2 ) + rj 



sin(ki z di)cos(k 2z d 2 ) + T] + cos(ki z di)sin(k 2z d 2 ) — 77 sin(k 2z d 2 ) 
sin(ki z di)cos(k 2z d 2 ) + rj + cos(ki z d\) sin(k 2z d 2 ) + rj~ 



-rj + cos(A:i,(ii)sin(K2 z ^2) — 1? s,in(k 2z d 2 ) 
- rj + cos(ki z d\) sin(k 2z d 2 ) + rj~ sin(k 2z d 2 ) ' 



(8) 



where the phase \j/ e = K e d and the K e is the effective propagation constant. The subscript 1 and 2 refer to the layer A 
and B, respectively. The 3?\ is the generalized impedance of the equivalent layer A. The z-component wave vector k\ z 
and k 2z are determined from the dispersion relations in Eq. (0. The d\ is twice the thickness of layer A; and the d 2 is the 
thickness of layer B. The period of the symmetric unit cell is d = d\ + d 2 . Other parameters in Eq. ([8]) are given by 



± _ 1 f £2xku £ixhz 
2 \£i x k 2z £2xk\ z 

± = I ( V2xh z L li\ x k 2z 
2 \Huk2z toxku 



for TM , 
forTE. 



(9) 



Although some MMs are mesoscopic media, nevertheless we restrict our consideration to those satisfying rf«l, after 
tedious derivation, Eq. (H)) can be simplified to 



TM: 



TE: 
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(10) 
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Here, the Ej and jUi (j = x,y,z) are the bulk values of the effective permittivities and permeabilities, respectively. They 
are related to the material parameters of the equivalent layers A and B through 



M P = 



d\E\ p + d 2 E 2p 

d ' 

d\jX\p + d 2 jX 2p 



£7 = 



Mz = 



d £\ z E 2z 

d\ £ 2z + d 2 £\ z 
dHi z H2z 



(11) 



d\\i 2z +d 2 \L\ z 



where p = x,y. Equation ( fTTT i is a simplification of the most general linear case treated by Lakhtakia lPTBI . In practice, 
the material parameters of the equivalent layers are unknown. One can only measure scattering parameters. In most 
experiments 3ft = <2o, combining Eqs. (O and ( [Tol l, we obtained the retrieval formula: 
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where 



X = e b n h sm 2 e, Y m = —Scos 2 9, Y e = —Scos 2 9, S~ 
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(13) 
(14) 



where 9 is the incidence angle. The Et, and jj,t, are, respectively, the background relative permittivity and permeability. 
The retrieval formulas in Eq. (TT2b provide four straight lines (Ym, Y m , Ye, Y e vs. X), two for each polarization. The Ym 
and Ye represent the dispersion lines for TM and TE polarizations, respectively. The Y„, and Y e are the corresponding 
impedance and admittance lines. These straight lines are easy to implement experimentally. After measuring the 
scattering parameters at several incidence angles and plotting the data according to Eqs. (fT2l-([T4li. use linear regression 
technique to calculate the slopes and Y-intercepts of the four lines. From the slopes and Y-intercepts, the six effective 
material parameters, Ej and ~jXj (j = x,y,z), can be retrieved simultaneously. Let Y M and Y% represent the Y-intercepts of 
the two lines in TM polarization, Sm and 5„, the corresponding slopes of the lines; whereas Y E \ F e °, Se, and S e are the 
corresponding quantities for TE polarization. Thus, 



TM: 
TE: 
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(15) 



where the ± sign in Eq. (fT3T > can be fixed by requiring the imaginary part of refractive index greater than zero, i.e. 
3(n m ) >0,3(«<;) > 0, and the real part of impedance and admittance greater than zero, i.e. 3i(Z) >0,and3i(Y) >0for 
passive media lfTTI . 



4. Resolving phase branch 

There are two issues regarding the branch determination in Eq. (fl4l , the ± sign in front of the inverse cosine and 
the phase branch m. The inverse cosine cos _I (-) represents the principal value: < cos~'(-) < %. Knowing only the 
cosine value, the phase angle cannot be determined. One need the information of sine, which can be obtained from the 
imaginary part of the cosine as following: 



I -sl 



-S 2 



25, 



: cos(p = cos(<p r + /<p,-) = cos (p r cosh (pi — i sin (p r sinh q>i , 



(16) 



where <p, and (pi are real. For passive medium, (pi > 0, and thus sinh 0, > 0. Therefore, the ± sign in front of the inverse 
cosine in Eq. ( fT~4T > can be resolved from the imaginary part of A: 



<P 



2%- 



if3(A)<0 
if3(A)>0 



(17) 



After determining the ± sign, next step is to resolve the correct phase branch m. This part can be very confusing when 
negative refractive index might be involved. Notice that among the four retrieval lines in Eq. (fT2l) . only the dispersion 



lines, Ym ~ X and Ye ~ X, are functions of the phase branch m, the impedance and admittance lines are independent of 
m. Employing this property, we provide three methods that might help to resolve the correct phase branch m. 

Method- 1, Using e x and ~p x . For the TM polarization in Eq. (fTZt . the e x can be obtained either from the intercepts 



Y° 



as £ x = -M- or from the slopes as e x = — — , and the two results should be close. When the branch is wrong, the two 



results can be significantly different. Similarly for the TE polarization, the \i x can be obtained from either fx x 
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— . Hence, the correct branch m is the one that minimizes the absolute value of the differences, i.e. 

Se 



TM: 



TE: 



mm 



mm 



Y°(m) S M {m) 



Y°(m) S E {m) 



Y? 



Mv 



£y 



Mv 



ra = 0,±l,±2, 



m = 0,±l,±2,- 



By 



(18) 



Method-2, Using -=- and =- : Alternatively, we can use ■=- for TM and z±- for TE as criteria to determine the correct 
e z H z E- n 



phase branch, i.e. 
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(19) 



Method-3, Using £ z pL y and £ y pL,: The third method is to use £-/X v for TM and EyPL. for TE as criteria to select the 
correct phase branch. Then, the algorithm becomes 
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TE: 
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(20) 



In above three methods, the branch number m predicted by the first two methods is always consistent for all the 
frequencies in the several examples we tested. The third method predicts the same result as the first two for most of the 
frequencies, but sometimes it can be different by ± 1 at the edge of phase transition frequencies in resonant regimes. Note 
that before applying Eqs. (fT8l . (fT9l , or d20l l to resolve the correct phase branch, the ± sign in front of the inverse cosine in 
Eq. (fT4b must be determined first. From the several examples we tested, using above branch-resolving techniques, there 
is no need to recourse adjacent frequencies in determining the correct phase branch. Recoursing adjacent frequencies 
can become confusing when both positive and negative refractive indices are present in the same frequency band. 



5. Discussion 

In this section, we provide an example to show how to implement the graphical method from the scattering parameters. 
Typically the background is lossless and thus, the X variable in Eq. (fT2l is real; whereas the slops and Y-intercepts are 
usually complex. The real and imaginary parts of the retrieval lines should be plotted separately as shown in Fig.|2]which 
was calculated from the scattering matrix. Drude model is used for the effective material parameters of the equivalent 
layers A and B, 
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n = i 
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(21) 



f 2 -f} r +iyf " f*-fZr+irf' 

where f ep = 30THz, f mp = 20THz, and y = 3THz. The e x and fi x are described by Eq. OTb with resonances at 
f er = 20THz and f mr = 25THz for the layer A and at f er = 35THz and f mr = 37THz for the layer B. The other 





(a) TM polarization: dispersion lines (left panels) and impedance 
lines (right panels). 



(b) TE polarization: dispersion lines (left panels) and admittance 
lines (right panels). 



Fig. 2. Retrieval lines at the frequency 30THz calculated from the scattering matrix of a 
single layer of unit cells at seven incidence angles (denoted as circles) uniformly distributed 
from to 30 degree. The horizontal axis is sin 6. Background medium is vacuum. Top 
panels: real part. Bottom panels: imaginary part. 



parameters for A: e y = e A — 0.3, £- = E x + 2, fX y = fX x — 0.5, and fi : = 1. The other parameters for B: E y = e x — 0.8, 
E z = £, — 0.5, H y = jJ. x + 0.2, and jj. z = ji x — 0.6. The thickness is 240nm for the layer A and 320nm for the layer B. 
Thus, the period of the unit cell (ABA) is 800 nm. Figure[3]shows the correct phase branch m predicted for each frequency 
in the regime of interest for the thickness of one (a) and six (b) periods. In our example since d <C A, for one unit-cell 
thickness, most of the frequencies are within the fundamental branch m = except for the regime of negative index of 
refraction [see Fig. 0a)] where m = — 1. For the six-period thickness, the phase branch jumps to m = 1 in the frequency 
range 18 ~ 20THz due to the high valves of the positive refractive index in this regime[see Fig.Ua)]. 

Shown in Fig. @] are the effective index of refraction, impedance, and admittance retrieved from one unit cell. These 
values are the same as those retrieved from the six unit cells. In this example, both TM and TE polarizations contain two 
frequency bands where the effective index of refraction is negative. As illustrated in Fig. [5J these two negative bands 
are accurately captured by the branch-resolving techniques provided in Sec. |4] The six retrieved effective permittivities 
and permeabilities are shown in Figs.|5]and|6]in the blue-solid curves, which are obtained from the graphical-retrieval 
and phase-unwrapping scheme, i.e. Eqs. (JT3J, ( TTTb , and d20l ). For comparison, we also show the effective permittivities 
and permeabilities calculated from Eq. (fTTT > in Figs.[5]and[6]in the green-circle curves. The retrieved effective material 
parameters agree very well with the results calculated from the effective medium theory. These values represent the bulk 
values of the material parameters since they are extracted from the symmetric unit cells. We have successfully applied 
the graphical method and phase-unwrapping techniques to our recent experiment ifTTl . From the experimental point of 
view, the straight-line graphical method is more accurate than the methods using single data point. The linear regression 
technique is based on collective data points measured at several incidence angles. It has an averaging effect and thus 
reduces the uncertainty of measurements. Using several incidence angles can also help to resolve the phase ambiguity. 
After determining the effective parameters of the bulk material, Eq. (fTTT) can be used to recovery the effective parameters 
of the equivalent layers if we are interested in. This extra benefit may help to improve the unit-cell design. 

When the tensors in Eq. (HJ rotate about the z-axis without rotating the coordinate system, the off-diagonal elements 
in the x-y plane will not be zero (e xy ^ and ji xy ^ 0), and thus the x and y components of the electromagnetic fields 
will be coupled. The dispersion relations in Eq. (O will no longer be valid. Without the proper modification, the current 
parameter retrieval scheme cannot be applied to this scenario. The original classification of the left-handed and right- 
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(a) One-period thickness. 
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(b) Six-period thickness. 



Fig. 3. Branch number predicted by the method-3 [Eq. d20H vs. frequency. Since rf<A, 
most of the frequencies are in the fundamental branch m = 0. The m = — 1 indicates the 
frequencies of the negative refractive index [n < 0, see Fig. |4(a)| . The m=\ corresponds 
to the frequencies of the high valves of the positive refractive index [see Fig. |4(a)| . Top: 
TM polarization. Bottom: TE polarization. 
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(a) Retrieved refractive index for TM polarization (left panels) and 
for TE polarization (right panels). 
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(b) Left panels: Retrieved impedance for TM polarization. Right 
panels: Retrieved admittance for TE polarization. 



Fig. 4. Effective refractive index, impedance, and admittance retrieved from a single unit- 
cell layer. Top panels: real part. Bottom panels: imaginary part. 



handed electromagnetic materials can cause confusion with chiral materials which are important class of electromagnetic 
materials OH- To avoid such confusion, we can use phase-velocity characteristics. When the phase velocity and the 
time-averaged Poynting vector form an acute angle, the medium is positive-phase-velocity characteristics. When the 
phase velocity and the time-averaged Poynting vector form an obtuse angle, the medium is negative-phase-velocity 
characteristics. 
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(a) The effective E x (left panels) and ~p x (right panels). 
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(b) The effective £,, (left panels) and ~p (right panels). 



Fig. 5. In-plane values of the retrieved material parameters. Upper panels: real parts. Lower 
panels: imaginary parts. 
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Fig. 6. The effective e- (right panels) and fl. (left panels). Upper panels: real parts. Lower 
panels: imaginary parts. 



6. Conclusions 

In conclusions, we have shown that symmetric unit cells can often be constructed in metamaterials regardless the number 
of sub-layers and the complexity of each sub-layers in the original unit cells. The graphical retrieval method and phase 
unwrapping techniques presented here might be a useful tool for metamaterial designs. 
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